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Abstract 



A recent renormalization group approach to a modified Bak-Sneppen model 
is discussed. We propose a self-consistency condition for the blocking scheme 
to be essential for a successful RG-method applied to self-organized critical- 
ity. A new method realizing the RG-approach to the Bak-Sneppen model is 
presented. It is based on the Monte-Carlo importance sampling idea. The 
new technique performs much faster than the original proposal. Using this 
technique we cross-check and improve previous results. 
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I. INTRODUCTION 



Biological evolution has been a candidate for self-organized criticality (SOC) PJ for a long 
time. The Bak-Sneppen (BS) model was the first model describing biological evolution 
as a self-organized critical phenomenon. A new renormalization group (RG) approach to 
calculate critical exponents of the BS-model has been presented first in |J. This note presents 
a Monte Carlo technique to realize the approach of [[J and provides thereby an easy way of 
cross-checking the novel RG-method. It accelerates and simplifies the calculation by using 
a Monte-Carlo inspired technique. We are able to improve the results given in ||. The 
paper is organized as follows: in section [TJ we redefine the BS-model and review the main 
ideas of the RG-approach. The main tool of this approach is the Run- Time-Statistics (RTS) 
technique, first presented in M. We discuss it in section [TII[ Section [TV] introduces our new 
technique to perform the RG-approach more efficiently The fixed point properties are given 
and basic critical exponents are calculated. 

II. THE RG-APPROACH TO THE BS-MODEL 

Following Ref. we modify the BS-model and define a left-or-right L/R-BS-model: 
Consider N real variables (t) = Of } ), % G {1,...,JV} with values 0.0 < < 1.0. At 
each time t determine index i(t) indicating the smallest value fa(t) among all fy. Replace 
the value fa^) by a new value x with probability a(x), J a(x)dx = 1. Choose with equal 
probability the left or right neighbor of the active site i(t) and replace it by a new value 
y with probability b(y), J b(y)dy =10. After equilibration almost all variables <f)j have a 
higher value than some p c . A p-avalanche is defined to start at time ti with the minimal 
value (pittA = Pi lasts as long as consecutive minimal values are smaller than p and stops 
at tf if 4>i(tt) > P- The temporal extension is given by s = tf — t^. The spatial extension 
£ is defined as the largest extent of active sites involved, £ = maxt^^^ \i{t\) — 2(^2) |- 
Avalanches keep themselves running by generating variables smaller than p. If this number 
of variables smaller than p, n t (p) = Ya @(/°~ 0i )> decreases with time during a p-avalanche 
the avalanche is called sub-critical and will die out. On the other hand, if n t (p) increases, 
the avalanche is called super-critical and will last forever. If n t (p) stays constant we observe 
a critical avalanche. Critical avalanches obey power laws for their spatial and temporal size, 
•P(t(£) ~ £~ Tf = and P(s) ~ s~ T respectively. The connection between spatial and temporal 
size is given by s ~ The critical exponents r^, r and z are connected via the scaling 
relation zt = — 1 + z. Extensive numerical studies have been performed to calculate these 
exponents for the original BS-model |6|-§. 

The goal of a RG-approach is threefold: calculate critical exponents in an independent 
analytical or semi-analytical way; prove the attractive nature of the critical state; demon- 
strate the concept of universality. The RG-approach considers small avalanches as objects 
to be integrated out and provides a mapping of larger avalanches onto smaller ones. This 
goes in parallel with a modification of the dynamical rules at coarser scales. 

To be more precise, let us denote the fine scale with index (/), the next coarser scale 
with (I + 1). Block variables are $/, fine variables are fa. Using a block factor of two, the 
simplest block transformation conserving the spatial avalanche structure is 
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<lf > = min^^} ■ (2.1) 

Then the dynamics on the coarser level is again based on the selection of the minimum. 
To perform one update on the coarse grid, i.e. updating 3>/(t) and either the left or right 
neighbor I(t) ± 1, we have to consider a process with four neighboring variables on the finer 
scale being updated. Thus a dynamical blocking of time occurs and we have different time 
scales on the fine and coarse level. The mapping connecting both time scales is dynamic, 
i.e. it depends on the stochasticity of the temporal evolution. This will provide us with 
a modified dynamics based on new distributions a" +1 '(x) and b" +1 ^(x). The initial and 
final state of one updating step on the coarse level is obtained by applying equation ( |2.1|) 
to the initial and final configuration of the corresponding fine level process. The spatial 
and temporal blocking of the dynamical process is schematically shown in Fig. |I[ The 
block transformation and the modified dynamics have to be designed such that it makes 
statistically no difference for the final state $^ if we apply first T fine updates and then 
the block transformation or if we apply first the block transformation and then perform one 
coarse update using the coarse dynamics: the blocking diagram has to be commutative for 
a successful RG-approach. 

Assume that we already know the probability distribution a« +1 \x) and ^ +1 )(x) of 
variables $i(t) an d &i(t)±i m the final state on the coarse level. It is then possible || 
to perform a variable transformation x — > y(x) using y(x) = Jq a" +1 ' (x')dx' leading to 

a ('+i)( y ) = Q(y)Q(l- y ) and b( l+1 \y) = fSfl . This reduces the possible RG- 

proliferation to the distribution b(x) . All other rules (search for the minimum, choose left or 
right neighbor, replace the value of the minimal variable with a uniformly distributed one) 
are invariant under the RG-transformation. 

Scale invariant behavior is expected only for avalanches that are critical, i.e. <t>il.) — pP ■ 
Therefore one has to consider only critical processes for the mapping of a sequence of fine 
updates onto one coarse update. Only processes a with four neighboring sites being updated 
after T a time steps are relevant. Processes that stop earlier because 4>f^ > pf for a > u 
do not contribute since they represent small fluctuations with avalanche size S < 3 and 
are neglected in the RG step. The parameter pp is determined from the condition that 
n t (p^\b^) should stay constant for critical avalanches, see eq. ( |Q| ) below. 

III. REALIZATION OF THE RG-APPROACH USING RUN-TIME-STATISTICS 

We now want to discuss the Run- Time-Statistics approach developed in Q. Denote 
the probability distribution of <pi at time t by pi tt (x). The conditional probability that the 
minimal variable is located at site i is then 



?4in,i0) = PiA x ) II / PjAy) d y ■ 



For the blocking procedure we are interested in the probability that (pf 1 is the minimal site 
and that the avalanche does not stop, i.e. </>f < p, which is 

Vi,t(p) = [ P Pmm,i( X ) dx ■ ( 3 - 2 ) 
JO 
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Exploiting the information about the position of the minimal value (pi, (pi < p, the probability 
distribution Pj t t(x) of all other sites j modifies in the next time step to 

Pj.t+i(x) = Npj. t (x) f p it (y)G(x - y) 
Jo 

X II / PkAz)dzdy . (3.3) 

jV is a normalization factor. The probability distribution of the minimal site i(t) and its 
left or right neighbor change according to the rules of the L/R-BS-model to 

Pi(t),t+i{x) = 9(x)6(l - x) (3.4a) 
Pi(t)±i,t+i{x) = K x ) ■ (3.4b) 

Note that equation ( |3.3| ) is different from the corresponding equations in PH . This leads 
to different results of the fixed point properties, see below. Using the iterative rules (|3.3|) 
and ( |3.4| ) we can calculate the weight v a of a given process a contributing to the blocking 
procedure 

■y To, 

V «(P) = 7^T[Via(t),t(p) ■ ( 3 - 5 ) 

Z " t=l 

The final probability distributions in the four updated variables contribute to the blocked 
distributions a and b at the next coarser level according to the weight of process a. Using 
the variable transformation described in [fj], a(x) may be rescaled to a uniform distribution 
in the interval [0, 1] and the next renormalization procedure can be iterated. 



IV. THE MONTE-CARLO REALIZATION OF THE RG-APPROACH 

To realize the renormalization group approach using the RTS-technique, one has to 
calculate the weights and final probability distributions of all contributing processes for 
each iteration during the computation of p c , using equation (|4.3|) below. Then one has to 
perform the same task for all paths using the value p c . Finally the block transformation ( [2.1|) 
has to be applied. However, there is a huge amount of contributing processes: up to a length 
of T a < 20 we estimated a total number of O(10 12 ) contributing processes. It is also clear 
that most of them, especially of the longer processes (e.g. the ones with the minimal site 
at the same place for a long time) are of very small probability. This observation reminds 
to the problems in "simple sampling" Monte Carlo algorithms, calculating properties and 
probabilities of all possible states in phase space. It is more convenient to use an " importance 
sampling" method generating the contributing processes a according to their weight v a . 
Then all generated processes of the ensemble contribute with equal probability. In other 
words, we evaluate RTS-integrals using the Monte-Carlo importance sampling method. 

To generate a single relevant process a. using distribution b^ l '(x), we start with a variable 
4>i a (Q) = p, and apply the rules of the L/R-BS-model until four neighboring sites are updated. 
For this process we count the number of variables smaller than p after the first time step as 
well as in the final state 
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n% = EMP ~ <Pt Ta) ) ■ (4-lb) 

If n[ a) = or rij^ = 0, the generated process is not a relevant process, since it stopped, 
i.e. > p, before four neighboring sites are updated. Using this notation we can easily 
write down n t =i(p,b^) and ^=t flnal (p, b">) obtained from the generation of an ensemble of 
iV processes [l(J 

n t=1 (p,& (0 ) = ^E^ a) (4.2a) 

a 

Since we are interested only in critical avalanches we have n t (p^ , b^) = const, and obtain 
p(y as solution of 

n 1 (pW,6W)=n tflnal (p« ) ^) . (4.3) 

Using the secant method it converges in about five iterations. Once we know the critical 
pV> at level / for given we are able to evaluate all critical relevant processes a starting 
with 4>i a (o) = p\P for their contribution to the blocked distributions d^ +1 \x) and W +1 ^(x). 
A relevant process a contributes with 

a« +1 \x) = 8(x - m a ), m a = min{4* a) , ^li) (4.4a) 
b { l +l \x) = 8(x - n a ), n a = mm{^\ • (4.4b) 

The first active site i a (t=o) has index 2j or 2j + 1 and the pair (2k, 2k + 1) is the right or left 
neighbor of (2j, 2j + 1). The four indices 2j, 2j + 1, 2k, 2k + 1 have been updated in process 
a. For an ensemble of N' relevant processes a, each occurring with a probability u a (p^), 
we obtain 

~a {l+1 \x) = ^j:~a^\x) (4.5a) 

a 

^w-kP'w • ( 4 - 5b ) 

J V Q 

To these preliminary distributions we apply the variable transformation described above, 
thereby switching back to a uniform distribution a^ +1 ^(x) and a transformed distribution 
b^ l+1 \x). Now the described procedure may be iterated. 

An ensemble of processes a, each occuring with probability u a , is generated by applying 
the rules of the L/R-BS-model to processes starting with 0» a (o) = p as long as not four 
neighboring sites have been updated. Once this happens, we stop the process and save 
its initial and final state. We reject non-relevant processes. Then the process a occurs 
automatically with its correct weight u a without the need of an explicit calculation of v a 
using RTS. Then we start the same procedure again for the next process that is generated 
completely independently from the previous one. 
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Being considerably simpler our method is able to reproduce the results of [|J with good 
accuracy in about one hour CPU-time on a workstation. Moreover, there is no need to 
extrapolate the results from T max = 20 to T max = oo as done in 0. In our approach we 
effectively consider T max = oo, since we allow arbitrary long relevant processes. 

On the other hand, since we generate a finite number TV of processes, we observe statis- 
tical fluctuations in the resulting distribution b^(x), whereas the method || gives in this 
sense an exact result (neglecting small errors due to the numerical integration routines). But 
this is a relatively small drawback compared to the advantage of speeding up the method, 
thereby having an effective method at hand to cross-check previous results. 

Starting with a distribution b^°\x) = 0(x)O(l — x) at the finest level / = 0, we observe 
that the distribution b^(x) converges very fast. It reaches its fixed point shape already 
after one renormalization step, see Fig. fj. Correspondingly the value of p^ converges also 
very fast to its fixed point value p* = 0.5954. For the calculation of b"'(x) we evaluated 10 9 
relevant processes at each renormalization step. For the determination of we evaluated 
10 7 processes per iteration of the secant method and renormalization step. The results of 
differ significantly which is probably due to the error in one of the RTS-iteration equations 

The probability distribution P*(T) for the length T of relevant processes, corresponding 
to the fixed point properties b* and p*, decays exponentially: P*(T) ~ e -constxT_ rj-^g 
minimal length of a relevant process is T = 3. At T = 30 the probability is of order 10~ 5 . 

To compute the critical exponents from the fixed point properties we follow the discussion 
of and use equations (5) and (6) given therein. For details the reader is referred to ||. 
From the probability distribution P*(T) we obtain the exponent z = 2.336(1). For the 
calculation of exponent r^, describing the spatial scaling of avalanche sizes, the probability 
K*, that the site with minimal value is always the same (until the avalanche stops), has to be 
known. The value of K* may be obtained again using our Monte-Carlo approach. We have 
to count the fraction of processes with the minimum at the same site at each time step until 
they stop due to the condition <pig(T ) > p\- Using the fixed point distribution b*(x), we find 
K* = 0.1827(1), from which follows = 1.2911(1). Using the scaling relation, connecting 
T^, r, and z, and our result for the dynamical exponent z we obtain for r = 1.1246(1). 

The values of exponents obtained via the renormalization group approach and via nu- 
merical simulation are compared in Table |. Due to the probabilistic nature of our approach 
the fixed point distributions fluctuates a little. The relative error in p*, z, K*, and r 
induced by this fluctuation has been estimated of order of 10~ 3 . The results obtained by 
numerical integration of the RTS-equations given in || differ significantly due to the error in 
one of the RTS-equations. The results of the Monte Carlo RG-approach show a coincidence 
up to 4% with numerical values for the exponents of the original BS-model. It is assumed 
that the L/R- version is in same universality class as the original BS-model. 

V. CONCLUSIONS 

We have introduced a new Monte Carlo renormalization group method for the L/R- 
Bak-Sneppen model. Previous results for basic critical exponents are improved. We have 
proposed a self-consistency condition for the blocking diagram in RG-methods applied to 
SOC systems. 



6 



ACKNOWLEDGMENTS 

I would like to thank G. Mack and Y. Xylander for stimulating discussions. This work 
was supported in parts by the Deutsche Forschungsgemeinschaft and the Studienstiftung 
des deutschen Volkes. 



7 



TABLES 



TABLE I. Comparison of results for critical exponents using the renormalization group ap- 
proach and simulations. The renormalization group approach applies to the L/R- version of the 
BS-model. RTS-approach results are given in |J 





MC approach 


RG method 

RTS approach 


Simulation 
original BS 


T 


1.1246(l) a 


1.1204 a 


1.08(1) 




1.2911(1) 


1.2766 




z 


2.336(1) 


2.2975 


2.43(l) b 



a These results have been obtained from the scaling relation. 
b This result is given in [fjffl. 
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FIG. 1. Concept of blocking in space and time using a block transformation BT. T updates 
on the fine level I are performed until four neighboring variables are changed. The result of this 
process defines a contribution to one time step on the coarse level / + 1. 
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FIG. 2. The fixed point distribution b*(x) is observed already after one renormalization step. 
For each RG-step 10 9 relevant processes have been generated. 
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